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STATEMENT REGARDING FEDERALLY SPONSORED 
RESEARCH OR DEVELOPMENT 

The present invention was supported in part by grant number ROl CA72023 from the 
National Institutes of Health (NIH/NCI). The U.S. Government has certain rights in the 
5 invention. 

FIELD OF THE INVENTION 

The present invention relates generally to medical imaging. More particularly, the present 
invention relates to methods for detecting and classifying a structure of interest in medical 
10 images to differentiate for instance (pre-)cancerous tissue from normal tissue. 

BACKGROUND 

Computed tomographic colonography (CTC) was first suggested in the early 1980s as a 
15 potential method for mass screening of colorectal cancer, the second leading cause of 
cancer deaths in the US. CTC was first realized in the 1990s following the rapid progress 
in computed tomography (CT) and in digital computing. CTC is a minimally invasive 
method that involves the steps of CT imaging the whole abdomen and pelvis after cleansing 
and air insufflation of the colon. Since the first realization, several studies have been 
20 conducted assessing the performance of CTC, mostly based on a radiologist's visual 
examination of either two-dimensional (2-D) CT images or three-dimensional (3-D) virtual 
colonoscopic views, or both. Most efforts have been directed toward developing better 
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visualization and navigation techniques, such as rendering, colon wall flattening, flight path 
planning algorithms, and user interface design. Recently some research has focused on 
developing computer-aided detection (CAD) methods for the identification of colonic 
polyps in 3-D CT data to improve the accuracy and efficiency of CTC. In these 

5 identification approaches, the 3-D geometrical features of polyps are extracted and used for 
their detection and identification. Mir et al. reviewed a set of methods proposed for shape 
description in CT images, e.g., moments, medial axis transforms, splines, curvature, Fourier 
descriptors, AR (Auto-Regressive) modeling, and statistical approaches (See A. H. Mir et 
al, "Description of shapes in CT images: The usefulness of time-series modeling 

10 techniques for identifying organs," IEEE Eng. Med. Biol. Mag., vol. 18, pp. 79-84, 
Jan./Feb. 1999). Summers et al. concluded that detection by shape analysis is feasible, 
especially for clinically important large polyps (See e.g. R. M. Summers etal, "Automated 
polyp detector for CT colonography: Feasibility study," Radiology, vol. 216, no. 1, pp. 
284-290, 2000; R. M. Summers et al, "Automated polyp detection at CT colonography: 

15 Feasibility assessment in a human population," Radiology, vol. 219, no. 1, pp. 51-59, 
2001). Paik et al. proposed to use a method based on overlapping surface normals to detect 
spherical surface patches along the colon wall that are likely to be parts of polyps (See e.g. 
D. S. Paik et al, "Computer-aided detection of polyps in CT colonography: Free response 
ROC evaluation of performance," Radiology, vol. 217(SS), p. 370, 2000; D. S. Paik et al., 

20 "Detection of polyps in CT colonography: A comparison of a computer aided detection 
algorithm to 3-D visualization methods," in Proc. 85 th Scientific Sessions Radiological 
Society of North America, vol. 213(P). Chicago, IL, 1999, p. 428). Yoshida et al. reported 
that geometric features extracted from small volumes of interest are effective in 
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differentiating polyps from folds and feces (See H. Yoshida et al, "Detection of colonic 
polyps in CT colonography based on geometric features " Radiology, vol. 217(SS), p. 582, 
2000), as well as characterizing colon wall surface geometry (See H. Yoshida et al, 
"Three-dimensional computer-aided diagnosis scheme for detection of colonic polyps" 
5 IEEE Trans. Med. Imag, vol. 20, pp. 1261-1274, Dec. 2001). Goktiirk et al fitted local 
spheres to the colon wall and based their detection on the existence of clusters of sphere 
centers (See S. B. Goktiirk et al, "A graph method for the conservative detection of polyps 
in the colon," in Proc. 2nd Int. Symp. Virtual Colonoscopy, Boston, MA, 2000). 

10 Most of these prior methods are rather sensitive (i.e. ability to detect positives), but need to 
be more specific (i.e. ability to detect negatives) for clinical applicability. The low 
specificity of some of the previously reported methods is generally due to the assumption 
that high curvature surface patches occur only on polyps. While it is true for instance that 
polyps have highly curved surfaces, so do some other structures, like haustral folds and 

15 retained stool. Radiologists reading these images use additional information to classify 
suspicious regions. For example, haustral folds are elongated structures, as opposed to 
polyps, which protrude locally from the colon wall. Stool may sometimes be identified by 
relatively inhomogeneous image intensity compared to polyps. However, if an automatic 
CAD method results in a low specificity manual examination of a (large) number of images 

20 corresponding to the CAD outputs is required to ensure proper detection. Such an 
examination is costly, time consuming and inefficient. Accordingly, there is a need to 
develop a method that would be capable of increasing specificity without sacrificing 
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sensitivity. Such a method could also be used to enhance and classify outputs of a high- 
sensitivity low-specificity CAD method to eliminate false positives only. 

SUMMARY OF THE INVENTION 

The present invention is a method capable of detecting and classifying a structure of 
interest with a high specificity without sacrificing the sensitivity of detection. The method 
could be used as a stand-alone method or as a post-processing method to enhance and 
classify outputs of a high- sensitivity low-specificity method to eliminate false positives. 
The method is based on representing changes in 2-D cross-sections of three-dimensional 
image data with a vector field, characterizing the topology of this vector field and using the 
characterized topology of the vector field for detection, identification or recognition tasks 
of a structure of interest. A structure of interest is a structure of interest to a radiologist 
such as polyps (e.g. colonic polyps), nodules (e.g. liver and lung nodules), lesions, or the 
like. 

More specifically, the method defines one or more image planes in a subvolume of a three- 
dimensional medical image. For each image plane, the edge displacement fields are 
computed for a plurality of slices that are defined over the axis perpendicular to its 
corresponding image plane. In general, one could define at least one image plane or a 
multiple set of image planes with each scrolling axis perpendicular to its 
respective/corresponding plane. In one exemplary embodiment, the image planes could be 
at least two mutually orthogonal planes. A more specific example could be the three 
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anatomical planes in which the scrolling axis is defined as the axis perpendicular to its 
respective plane, i.e. axial, sagittal, coronal. The edge displacement fields for all slices in 
an image field are combined to create an edge displacement field for that image plane. This 
combined edge displacement field is used to determine parameters. These parameters are 

5 then used to classify the structure of interest and determine whether the structure of interest 
is e.g. a polyp or a non-polyp. The classification could be based on the parameters from a 
single image plane or a combination of image planes. For instance, the three anatomical 
planes could be used to determine the parameters. In case three parameters are computed 
per image plane, one will end up with a 3x3 feature vector. The parameters could be 

10 further processed by taking for instance the average over the image planes for each 
parameter or by applying other functions or processes to combine the parameters and 
classify the structure of interest. 

BRIEF DESCRIPTION OF THE FIGURES 

15 The objectives and advantages of the present invention will be understood by reading the 
following detailed description in conjunction with the drawings, in which: 
FIG. 1 shows an example of generating an axial edge displacement field demonstrated 

using a sphere according to the present invention; 
FIG- 2 shows an example of three sequential axial images (smoothed for visual purposes) 
20 around a pre-detected structure and the associated EDF v w (*,)>) ; and 

FIG- 3 shows an example of axial EDFs of four cases: (a) positive (polyp) with [a; (3; d] = 
[-0.96; 0.29; 4.99]; (b) positive with [a; P; d] = [-0.99; 0.12; 5.15]; (c) negative 
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(non-polyp) [a; (3; d] = [-0.77; 0.65; 1.28]; (d) negative with [a; (3; d] = [-0.71; 
0.70; 0.72]. 



DETAILED DESCRIPTION OF THE INVENTION 

5 Although the following detailed description contains many specifics for the purposes of 
illustration, anyone of ordinary skill in the art will readily appreciate that many variations 
and alterations to the following exemplary details are within the scope of the invention. 
Accordingly, the following preferred embodiment of the invention is set forth without any 
loss of generality to, and without imposing limitations upon, the claimed invention. 

10 

The method of the present invention could be referred to as an "Edge Displacement Field 
Method" (EDF), a "Gradient Field Method" (GF) or an "Optical Flow Field Method" 
(OFF) applied to detection and classification of structures of interests in medical images 
(See for an overview of optical flow computation S. S. Beauchemin et al, "The 

15 computation of optical flow," Comput. Surv., vol. 27, no. 3, pp. 433-467, 1995). The 
present method is based on representing changes in 2-D cross-sections of three-dimensional 
(3-D) image data (e.g. axial, coronal and sagittal gray scale CTC data) with a vector field, 
characterizing the topology of this vector field and using the characterized topology of the 
vector field for detection, identification or recognition tasks of a structure of interest. The 

20 present invention could be utilized as a post-processing method that would refine the results 
of a high sensitivity, low specificity pre-detection by increasing specificity without 
sacrificing sensitivity. The present invention is, however, not limited to its application as a 
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post-processing method since it could also be used as an independent method without pre- 
identification. 



The medical images of the present invention are digital or computerized images such as, for 
5 instance, but not limited to, a CT, an MRI, a digitized X-ray, or any other medical image 
application that could be converted or rendered to a digital image. The medical images 
could be 2-D images used to construct a 3-D image or a 3-D volumetric image. A structure 
of interest is, for instance, but not limited to, a structure that contains pre-cancerous tissue 
or cancerous tissue. Examples of structures of interest to a radiologist are, for instance, 
10 polyps (such as colonic polyps), nodules (such as liver and lung nodules), lesions, or the 
like. However, even though the present invention is described with respect to medical 
images, a person of average skill in the art will readily appreciate that the present invention 
could easily be applied in any type of application where it is necessary to characterize a 
structure and determine whether this structure corresponds to a structure of interest. 

15 

The first step of the method is the EDF computation to represent the changes in the location 
of edges in the images (e.g. tissue/air boundaries) as one scrolls through the 3-D data. For 
instance, as shown in FIG. 1, let xy be an image plane perpendicular to the scrolling axis Z 
and determine the EDF. In general, one could define at least one image plane or a multiple 
20 set of image planes with each scrolling axis perpendicular to its respective/corresponding 
plane. In one exemplary embodiment, the image planes could be at least two mutually 
orthogonal planes. A more specific example could be the three anatomical planes and the 
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scrolling axis defined as the axis perpendicular to its respective plane, i.e. axial, sagittal, 
coronal. This specific exemplary embodiment of three mutually orthogonal image planes is 
further described infra. 



5 The EDF equation for the xy plane is: 

V/^.v z= , + ^ = 0 (l) 

10 where v z=k (x,y) is the EDF defined on the plane that is perpendicular to the z axis and is 
located at that z = k. I z=k (x,y) is the associated image, i.e. the attenuation coefficient 
function on the same plane. v z=k (x,y) represents the dislocation of the edge at (x, y) along 
the local gradient from z = k to z = k+l. v z=k (x,y) is computed for all k, i.e. for all slices, 
within the subvolume except at the boundaries. In one example, k could be 25, but k is not 
15 limited to 25 and could also be a higher or a lower number. In general, the number for k 
depends on the maximum polyp size of interest and the data resolution. The positive z 
direction was defined to be outwards from the center slice. This consistency is required as 
v z=k (x,y)for all k are summed and (optionally) smoothed to get a composite EDF, 
V z (x,y), associated with the current subvolume and the scrolling axis Z (see FIG. 1; 
20 v z (*,)>) can be considered as a 2-D representation of a 3-D subvolume of the xy plane 
scrolled over the Z axis). Thus it is assured that the edges of polyp-like structures move 
inwards on the plane perpendicular to the scrolling axis. The composite EDF could also be 
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(optionally) smoothed. An example of a smoothing kernel is a Gaussian (a = 2 mm) whose 
size could be limited to 2a. The steps could be completed for all three orthogonal axes (Z = 
[Axial, Sagittal, Coronal]) resulting in three EDFs that encode information regarding the 
structure of interest. 

5 

The next step is to characterize the vector fields for a single image plane or a combination 
of image planes. The characterization is based on determining some landmarks in the 
EDFs. For example, to characterize a single EDF, selected from v axial (x,y), 

Veronal * sagittal (*. , °™ P"«* n ° de ( PN > * V^ty, e.g. eight child nodes 

10 (CNs) could be determined. A PN is defined to be the minimum divergence pixel location 
on the EDF (the PN is marked with a square in FIGS. 1-3). In case a pre-detection method 
is used to detect a structure of interest the PN could be defined to be the minimum 
divergence pixel location in, for instance, a 4 x 4 mm neighborhood of the pre-detected 
point on the EDF. CNs are defined to be the pixel locations that are a distance (e.g. 4 mm) 

1 5 away from the PN on the streamlines incoming to the eight immediate neighbors of the PN 
(CNs are marked with a circle in FIGS. 1-3). FIG. 2 shows an example with three 
associated axial images and the associated EDF (*,)>). The parent node (PN) is 
marked with a square and the children nodes (CNs) are marked with small circles. It is 
noted that two of the eight CNs coincide with two other CNs. 



20 



Each of the EDFs vV ( U,y), v coronal (x,y) and v sagittal {x,y) could be used to compute 
parameters to classify the structure of interest. As indicated supra the classification could 
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be done using parameters computed from one image plane or a using the parameters 
computed from a combination of image planes. In one example two parameters, a and p\ 
are computed using the Jacobian matrix of the EDF at the PN (See e.g. J. Helman et al., 
"Representation and display of vector field topology in fluid flow data sets," IEEE 
Computer, vol. 22, pp. 27-36, Aug. 1989; or Y. Lavin et al., "Feature comparisons of 
vector fields using earth mover's distance," in Proc. Visualization '98, pp. 103-109, 524). 
The computation of those parameters is as follows: 



10 



J = 







dx 


dy 






dx 


dy 



a = 



3v. dv, 



dx dy 



(2) 



(3) 



P = 



a 2 -4 J 



a -4 J 



(4) 



15 Note that a and P carry information about the eigenvalues of the Jacobian matrix J . In fact, 
the characteristic equation of J is: 
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A 2 -aA+|J| = 0 (5 ) 



Furthermore, a is equal to the divergence of the EDF at PN. p is also computed at the PN 

and describes the circulatory behavior of EDF at the PN. The ratio of a to P uniquely 
defines the topology of a linear vector field at the PN so the normalized a and (3 



(normalized by yCC 2 + fi 2 ) are used as suggested by Lavin et al. (Y. Lavin et al., 
"Feature comparisons of vector fields using earth mover's distance," in Proc. 
Visualization'98, pp. 103-109, 524). 

Additionally, one could characterize the behavior of the incoming streamlines around the 
PN using the parameter d, defined as: 



e^(CN k CNj)e [0,7c] (6) 



where 0,/s are computed with respect to the PN's. d describes the spread of CNs around 
the PN. Thus, d is used to characterize the spread of CNs around the PN. 

FIG. 3 shows an example of four axial EDFs computed for two positive (polyp) and two 
negative (nonpolyp) cases to provide an understanding of the meaning of EDF 
characterization parameters visually. In agreement with the intuition, the EDFs 
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corresponding to positive cases have a PN with an a close to 1 (negative divergence) and 
P close to zero (small circulatory behavior), and a large d (streamlines well spread around 
the PN) 5 i.e., a star-shaped topology, unlike the EDFs of negative cases. 

Referring back to the exemplary embodiment related to three mutually orthogonal image 
planes, each parameter could be determined for the axial, coronal, and sagittal EDFs, 
resulting in a nine-dimensional feature vector for each subvolume considered. In one 
aspect, one could select to use the mean values of each parameter over three scrolling axes 

as the final feature vector / = [CC 9 P 9 d] where stands for averaging over axial, coronal 
and sagittal parameters. However, in another aspect one could use other functions and/or 
criteria of the EDF parameters instead of taking the mean of EDF parameters. 

The parameters or the feature vector could now be used to classify the structure of interest. 
In one example, a binary classification could be established to classify between a polyp 
versus a non-polyp using a Mahalanobis distance based linear classifier [See PC 
Mahalanobis, On the generalized distance in statistics, Proc. Natl Institute of Science of 
India 12:49-55, 1936). The Mahalanobis distance of a vector f to the mean vector m r of a 
population Tis defined as: 




(7) 
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where C r is the covariance matrix of T. This distance is a standardized measure that: 1) 
automatically accounts for scaling; 2) takes care of correlations between features; and 3) 
can provide linear and curved decision surfaces. For classification purposes, T represents 
the training set and f represents a sample from the test set £2. Referring to the subset of 
5 polyps in T as T l5 and the subset of nonpolyps as Tq, the binary classifier is defined as 
follows: 

r fm -r fm +&<0=>feQ 1 

£ m r, f > m r 0 1 

otherwise f €= (8) 

10 

where Q { and £i 0 refer to the subsets of polyps and nonpolyps in £2. An alternative to 
the Mahalanobis distance-based classifier is the use of SVMs, which minimize training 
classification error as well as generalization error. In general, the present invention is 
however, not limited to the selection of a particular classifier. 

15 

The present invention has now been described in accordance with several exemplary 
embodiments, which are intended to be illustrative in all aspects, rather than restrictive. 
Thus, the present invention is capable of many variations in detailed implementation, which 
may be derived from the description contained herein by a person of ordinary skill in the 
20 art. For instance the present invention has been described in the context of a method, those 
skilled in the art will appreciate that the method of the present invention is capable of being 
distributed in the form of a computer readable medium of instructions in a variety of forms, 
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and that the present invention applies equally regardless of the particular type of signal 
bearing medium used to actually carry out the distribution. In other words, the present 
invention is also a program storage device accessible by a computer, tangible embodying a 
program of instructions or means executable by the computer to perform method steps for 

5 detection and classification of a shape in a medical image as described supra. Examples of 
computer readable media include recordable type media such as disks and CD-ROMS and 
transmission type media such as digital and analog communication links. In addition, the 
present invention could be implemented and coded in different programming languages 
and/or packages. All such variations are considered to be within the scope and spirit of the 

10 present invention as defined by the following claims and their legal equivalents. 
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